arXiv:1501.01413v3 [astro-ph.CO] 15 Aug 2015 
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It is proposed to use exact, cosmologically relevant solutions to Einstein’s equations to accurately 
quantify the precision of ray tracing techniques through Newtonian N-body simulations. As an 
initial example of snch a stndy, the recipe in (Green & Wald, 2012) for going between N-body 
results and a perturbed FLRW metric in the Newtonian gauge is used to study light propagation 
through quasi-spherical Szekeres models. The study is conducted by deriving a set of ODEs giving 
an expression for the angular diameter distance in the Newtonian gauge metric. The accuracy of 
the results obtained from the ODEs is estimated by using the ODEs to determine the distance- 
redshift relation in mock N-body data based on quasi-spherical Szekeres models. The results are 
then compared to the exact relations. From this comparison it is seen that the obtained ODEs can 
accurately reproduce the distance-redshift relation along both radial and non-radial geodesics in 
spherically symmetric models. The reproduction of geodesics in non-symmetric Szekeres models is 
slightly less accurate, but still good. These results indicate that the employment of perturbed FLRW 
metrics for standard ray tracing techniques yields fairly accurate results, at least regarding distance- 
redshift relations. It is possible though, that this conclusion will be rendered invalid if other typical 
ray tracing approximations are included and if light is allowed to travel through several structures 
instead of just one. 

PACS numbers: 98.80.-k, 98.80.Jk, 98.80.Es 

I. INTRODUCTION 

Supernovae observations are typically interpreted as 
indicating an accelerated expansion of the Universe 
[1, 2]. Such interpretations, leading to the inclusion of 
a cosmological constant into the cosmological standard 
model, are predominantly based on redshift-distance 
relations valid only within the Friedmann-Lemaitre- 
Robertson-Walker (FLRW) models [3]. Since the 
Universe is not exactly homogeneous and isotropic, it 
is important to know how deviations from an exact 
FLRW universe affect light propagation. By using 
exact inhomogeneous solutions to Einstein’s equation 
it has e.g been shown that inhomogeneities can affect 
light propagation such that it may look as though the 
Universe is undergoing an accelerated expansion even 
though it is not (see e.g. [4-6]). As of yet, no convincing 
results have been obtained indicating that this is in 
fact what happens in the real universe though. On the 
contrary actually, studies indicate that randomization 
and statistical averaging diminishes any deviations from 
FLRW results [7, 8] 

Even if inhomogeneities cannot explain the seeming 
accelerated expansion of the Universe, being able to 
quantify even small effects may be important for param¬ 
eter determinations based on observations made in an 
era of precision cosmology. 
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^ Note added after publication: it has later been pointed out to 
us, that though this seems to be the case in general when cosmic 
backreaction is vanishing (see e.g. [35, 36]), it is not necessarily 
the case when the models studied have non-vanishing backreac¬ 
tion (see e.g. [37]). 


Newtonian N-body simulations are an important 
tool for studying structure formation and are widely 
accepted as giving a correct description of the structure 
formation going on in the real universe. Comparing 
results from N-body simulations with observations is 
thus essential for using survey data to put restrictions 
on cosmological parameters determining the standard 
model of cosmology. Unfortunately, Newtonian mechan¬ 
ics cannot be used to study exact general relativistic 
light propagation since this requires knowledge of the 
exact metric. Typically, light propagation through an 
N-body simulation is thus studied through ray tracing 
techniques (see e.g. [9-12]), where bundles of light are 
propagated through a simulation box via the background 
metric. If at all, only at a discrete number of lens planes 
is the light deflected based on linearized relativistic 
gravity. Several improvements of the basic ray tracing 
scheme have been developed and the accuracy of the 
methods have been studied through different approxima¬ 
tion schemes (see e.g. [12-15]). For example, in [15] the 
standard ray-bundle method (see e.g. [16]) was sought 
improved by using the first order metric, rather than the 
background metric, to describe null-geodesics. 

In an era of precision cosmology, it is pertinent to 
know exactly to what extent theoretical predictions of 
light propagation through an N-body simulation can 
be trusted and the work presented here is another step 
towards that goal. 

In [17] it is argued that the perturbed metric in 
the longitudinal gauge gives a good description of the 
metric of the Universe. A dictionary giving a recipe for 
going between Newtonian N-body results and this metric 
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is also given. The work presented here is based on that 
recipe. Several versions of the recipe are given in [17], 
varying from the most simple form to a very involved 
’’Oxford” version. The version which in [17] is denoted 
the abridged version is the one which is used here 
because it is that version which is most compatible with 
standard ray tracing techniques - and it is the simplest 
of the recipe versions to use with N-body simulation 
data. The metric is then simply the usual first order 
metric in the Newtonian gauge with the potential and 
velocity fields being those obtained directly from an 
N-body simulation. 

The compatibility between a first order perturbed 
metric and the metric of the real universe is not clear. 
It is advocated in e.g. [17] that it can be used, while 
others show more skepticism and e.g. point out dangers 
with using perturbation theory to justify itself (see 
e.g [18]). Since the perturbed FLRW metric in the 
Newtonian gauge combined with versions of the recipe 
of [17] are typically the basis of ray tracing through 
N-body simulations, it is important to study how 
accurate the use of this recipe really is. In order to do 
this, the metric and dictionary described above will be 
used to derive a set of ODEs which can be solved to 
obtain distance-redshift relations in the corresponding 
spacetime. To study the accuracy of the obtained 
relation, the ODEs are used on mock N-body data 
constructed from quasi-spherical Szekeres models [19] 
under the standard assumption that Newtonian N-body 
simulations accurately reproduce relativistic structure 
formation. Since the quasi-spherical Szekeres models 
are exact solutions to Einstein’s field equations, the 
exact distance-redshift relations of these models can be 
obtained and compared with those obtained by using 
the recipe. 

There are several reasons for choosing the quasi- 
spherical Szekeres models for the accuracy checking of 
the ODEs. First of all, the quasi-spherical Szekeres 
models are counted amongst the most realistic cosmo- 
logically relevant exact known solutions to Einstein’s 
equations. Further more, the spherically symmet¬ 
ric limit of the quasi-spherical Szekeres models, the 
Lemaitre-Tolman-Bondi (LTB) models, are known to 
be reproducible by Newtonian N-body simulations 
[20]. Thus, checking the ODEs with LTB models gives 
both a check of reproducing light propagation through 
exact non-linear relativistically evolving structures and 
through structures evolving according to Newtonian N- 
body simulations. It is not surprising that LTB models 
are reproducible by Newtonian mechanics since they are 
spherically symmetric and their structure formation is 
scale invariant. It is to obtain slightly more generality 
that also the non-symmetric quasi-spherical Szekeres 
models are used in the comparison. The non-symmetric 
Szekeres models show significantly more growth of 
structure than LTB models [21-24] and whether or not 
these can be reproduced by N-body simulations is still 
unknown but is the subject of ongoing work of one of 


the current authors. 

It should be noted that the perturbed FLRW metric 
is diagonal while the Szekeres metric is non-diagonal in 
spherical coordinates. This immediately implies a short¬ 
coming of the Newtonian gauge metric. However, there 
is no reason to expect the local metric of cluster-void 
structures in the real universe to be diagonal. Hence, the 
non-diagonal metric of the Szekeres model may actually 
be considered a strength for the purpose of the study. 


II. THE ANGULAR DIAMETER DISTANCE IN 
THE NEWTONIAN GAUGE 

In this section, the ODEs needed to obtain the angular 
diameter distance Da{z) in the Newtonian gauge are 
given. The ODEs are obtained following the procedure 
presented in [25]. 

The onset of deriving the sought equations is the 
metric. Since Szekeres models are dust models, the 
anisotropic stress vanishes and the perturbed FLRW 
metric in the Newtonian gauge can be written using a 
single perturbation field ip: 

da^ = —c^(l -I- 2ip)dt^ + a^(l — 2ip) {dr"^ + r^dd"^ + sm^{0)d(p‘^) 

= -Tdf + Rdr^ + Fde^ + Pd(p^ 
( 1 ) 

The metric is assumed to be exactly described in 
this form and will be referred to as the Newtonian 
gauge metric. It should be stressed, that though the 
metric above looks like a perturbed FLRW metric in 
the Newtonian gauge, there is a small difference: the 
above metric is an attempt to obtain a relativistic 
description of the ’’underlying” relativistic spacetime 
corresponding to a Newtonian N-body simulation. As 
such, according to the recipe of [17], the potential ip is to 
be obtained from the exact, non-linear density field from 
the corresponding N-body simulation, and not through 
perturbation theory. This point will be revisited in 
section HI C. 

As implied in the expression for the line element 
above, spherical coordinates have been used for this 
work. By writing the metric functions as T, R, F, P, 
the equations derived below are written in a generic 
form so that they are valid with any coordinate choice 
and actually any spacetime as long as the metric is 
diagonal. In appendix A, the equations are written in 
their full length in spherical coordinates using ip instead 
of T, R, F, P. 

Light moves along null-geodesics, so the geodesic 
equations are needed to describe the light paths. Letting 
a dot denote differentiation with respect to the affine 
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parameter A, the geodesic equations can be written as: 
-2Tk^ = 2fk*-Tt{k*f+ Rt{k'^f+ Ft{k^f 

( 2 ) 

2R¥ = -2Rk'^-Tr{k^f + Rr{k'^f + Fr{k^f + Pr{k'^f 

( 3 ) 


- 2Fk%k^^ - 2k%F^p^kP + F^pk%)- 

2Fk% - Te^ik^ - 2Tek^k% + R^kl^ + 2i?.efc"F„ + 

FMk'^f + ‘2Ffik<^k^^ + Pfic.{k‘^f + 2P^ek^kt 

( 10 )’ 


2FP = -2Fk‘^-Te{k*f + Reik’'f + F0{k^f + Pe{k'^f 

( 4 ) 

2Pk^ = -2Pk^-T4k*f+R^{k^f+F^{k^)^+P^{k'^f 

( 5 ) 

Subscripted commas followed by coordinates denote 
partial derivatives with respect to those coordinates. 

Considering a light bundle with infinitesimal cross 
section dS and solid angle element 6il, an ODE for Da 
can be obtained as follows (see e.g. [26]): 


6S = D\5D. 2d\-aDA = d\nSS = 20dX = k°‘^dX 

( 6 ) 

9 denotes the optical expansion scalar (see e.g. [26]) and 
a subscripted semi-colon followed by a coordinate denotes 
covariant differentiation with respect to that coordinate. 
Summation over repeated indexes is implied, with Greek 
indexes running over 0 — 3 and Latin indexes running over 
1 - 3. 

Inserting the appropriate Christoffel symbols (see ap¬ 
pendix A) into this equation, the following ODE is ob¬ 
tained: 


ADa 

dX 


1 

2 


DA{k^, 


t 

T 


R 

R 


F 

F 




( 7 ) 


Once this expression has been used to obtain Da along 
a geodesic, the luminosity distance can be found as 
Di^ = (1 -|- z)‘^Da and from this the apparent magni¬ 
tude etc. can be obtained. 

From equation (7) it is apparent that /s']., and 

k't are needed and these are obtained by solving the 
ODEs that appear when taking the partial derivatives 
of the geodesic equations: 


2P = -2P,J^ - 2P4fc(l - 2k^{P^,sk^ + p,pk^j 
-2Pk% - T^cik^ - 2r,^/c*fc[„ + i?.,>a(F)2 + 2i?,^fc"fc(„+ 
F^^^iky + 2F^fc®/c®„ + P,yky + 2Py^kt 

( 11 ) 


In order to ensure that the geodesics described by 
the above equations are null-geodesics, the null condi¬ 
tion and its partial derivatives shown below can be used 
to set the initial conditions: 


= -T{ky + R{ky + Fik^)"^ + p{ky = o ( 12 ) 

-ryef - 2TA:‘fc*„ + Ryky + 2i?Ffc(„ + 

Fyky -k 2FA:®fc®„ -k pyky + 2Pk'^k*^ = 0 


The objective of this work is to asses the accuracy 
of the equations presented above when applied to 
N-body data. In order to do this, the exact relativistic 
spacetime corresponding to the N-body data is needed. 
The study is thus conducted by constructing mock 
N-body data corresponding to a quasi-spherical Szekeres 
model and its underlying LTB model. The above equa¬ 
tions are solved with the potentials corresponding to 
these two ” data sets”. The results are then compared to 
those obtained by solving the equivalent ODEs derived 
from the actual metrics of these models. 


= 2T„fc* + 2TkyP^ + 2k\T^^pk^ + Ty%)+ 

2fc‘a - Tyky - 2Tye^ + Ryky + 2Ry^k^^+ 
F^ky + 2FAk<^k\ + pyky + 2Py^kt 

( 8 ) 

27? = -2A - 2Rky% - 2F(i?,„;3fc^ + R^^)- 

27 ?F„ - TrUky - 2.Tyk*^^ + R,yky + 2 i?,,FF„+ 
F^yky + 2Fy^k\ + p,yky + 2py^ki 

( 9 ) ’ 


III. QUASI-SPHERICAL SZEKERES MODELS 

In this section, the quasi-spherical Szekeres models 
are introduced and the ODEs needed to obtain Da{z) 
in these models are given. At the end of the section, the 
relation between the Szekeres models and the Newtonian 
gauge metric is discussed. 

The Szekeres models [19] is a family of exact, in¬ 
homogeneous dust solutions to the Einstein equations 
which in their general form have no killing vectors [27]. 
A typical coordinate system to describe the Szekeres 
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metric in, is the (t, r,p, q)-systeni related to the spherical 
coordinate system by a stereographic projection [28] 

p- P = S cot( 6 i/ 2 ) cos((/)) 
q — Q = S cot(0/2) sin((/)) 

The functions P, Q and S are defined below. 

The {t,r,p,q) coordinate system is particularly useful 
since it renders the metric diagonal: 


ds^ = 


— C dt + ^ - - ~ dV 


e — k{r) 

E{r,p,qy 


{dp^ + dq^) 


(15) 


In a general Szekeres model, E is given as 
E = where 

S', P and Q are continuous but otherwise arbitrary 
functions of r, and e S {—1,0,1}. The quasi-spherical 
Szekeres models with e = 1 reduce to LTB models when 
P, Q and S are constant functions. Only the quasi- 
spherical Szekeres models are considered in the following. 


Inserting the metric corresponding to the line ele¬ 
ment (15) into Einstein’s equation for a dust-filled 
universe containing a cosmological constant leads to the 
following two equations: 




A 


3c2 


(16) 


2M^r - 6M%^ 
c2pA^A,r-A^y 


p = SttG/c^ 


(17) 


The function M = M{r) appearing in these equations is 
a temporal integration constant depending on the radial 
coordinate and corresponds to the effective gravitational 
mass at comoving radial coordinate r. 


A. Model setup 

The procedure that has been used for obtaining spe¬ 
cific models follows that introduced by K. Bolejko and 


2 For notational convenience, the coordinates t,r,d and 4> of fho 
Szekeres metric will notationally not be distinguished from the 
coordinates of the Newtonian gauge metric. The coordinates of 
these two spacetimes are not the same though and mappings 
between the coordinate systems are discussed in section III C. 


described in e.g [21]. This method starts by specifying 
an LTB model with a density distribution that is later al¬ 
tered through appropriate choices of the dipole functions 
P, S and Q to create non-symmetric models. Below, a 
vanishing cosmological constant and k{r) < 0 is assumed. 

Equation (16) can be written as an integral equa¬ 
tion: 


r s dA 

ct-ctb{r)= / = (18) 

do y2M/A-k 

Here, only models with a constant time of the big bang 
will be used, and the constant value will be set equal 
to zero, i.e. tb(r) = 0. Introducing a parameter 77 , the 
solution to the integral equation can then be written in 
parametric form as: 

M M 

A = -^(cosh( 77 ) - 1 ), t = (sinh( 77 ) - (19) 


In the following, the term ’’background model” is used 
for denoting the FLRW model that the LTB model 
tends to at r —>■ 00 . In the models studied here, the 
background model is the Einstein de-Sitter (EdS) model. 


The coordinate covariance in r is eliminated by 
setting A{tis,r) = f": where Us is the time of last 
scattering in the background model determined by a 
background redshift z oi z = 1100 . 

To specify the LTB model completely, only one more 
function is needed. Here, that function is chosen to 
be M{r). In order to specify M, the function is split 
into two parts by writing M = 5M{r) -\- Mq. The 
latter term corresponds to a background term such 
that Mq = "*^^ 2 ° — , where is the time dependent 

density parameter of the background model with density 
PEdS- dM is determined by choosing the desired initial 
density field. This is done by splitting the density into 
a background part and a ’’perturbation” and writing 
equation (17) in the LTB limit as: 


PEdS + 5p = 


2((Mo),.+(5M,,) 

cAA^A^r 


( 20 ) 


Since M is time-independent, this equation can be 
integrated at t = tis to obtain M once Sp{tis,r) has 
been specified. When this is done, equation (19) can be 
solved for k{r) at t = Us and afterwards it can be solved 
at any time t to obtain A(t, r). 


The models used here are specified by 5p{r,tis) = 

— 10~^0.5PpEds{tis)oie^P°^^^^E^ with a = 5.5, cr = 0.3 
and P as defined in equation (17). This corresponds 
to a void with an approximate present time diameter 
of 160Mpc. The present time void profile in physical 
coordinates can be fitted well to the formula given in 
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Density at comovuig z = 0.035 



FIG. 1. Present time 2D density profile of the LTB model. 
The comoving coordinates are normalized at t = tu in units 
of 0.1 Mpc. 


Density at coinoving z = 0.035 



FIG. 2. Present time 2D density profile of the Szekeres model. 
The comoving coordinates are normalized at f = tis in units 
of O.lMpc. 


B. Geodesic equations and the Angular diameter 
distance formula in spherical coordinates 

A complete set of ODEs appropriate for studying 
redshift-distance relations in Szekeres models was de¬ 
rived in [25]. The coordinate system used there was the 
{t,r,p,q) coordinate system. For the purpose of this 
work it was found convenient to follow the procedure 
presented in [25] and re-derive the equations for quasi- 
spherical Szekeres models in spherical coordinates. 

The line element of the quasi-spherical Szekeres 
model in spherical coordinates is: 

(1 - k) 

^(S'2 cot2(^) -h 2S^r cot(^)(Q,^sin((^)-h 

Pr cos((/))) -I- P^ -I- Q\)]dr'^ + cot(^)[(5 r cos((/))- 

P_r sm{(j))]drd(j) — 2 —[Q^^. sin((/)) + P^ cos{(j))+ 
E ' 

S^r cot{-)]drd9 + A^dO"^ + siT?{9)d4>^ 

^ ( 21 ) 

The line element written in this form can also be found in 
e.g [21]. For convenience, a simplifying notation for the 
metric functions is used in the following and the ODEs 
will be given in a notation corresponding to the line ele¬ 
ment written as: 

ds^ = —dt^c^ + R{t, r, 9, (j))dr‘^ + 2^{t, r, 9, (j))drd(j)+ 

2Q(t, r, 9, (j))drd9 + F{t, r)d9^ + P{t, r, 9)d(j)^ 

( 22 ) 

The definition of the metric functions P, d), 0,P, P 
is seen by comparing equations (21) and (22). In 
spherical coordinates E is given by P = 2 sinpe/ 2 ) 

cos(0)+sin(0)[P r cos(^)+Q^t' sin(^)] 

“E" “ ^ 'S ^ ■ 


[29] which describes spherically averaged void profiles 
according to Newtonian N-body simulations (based on a 
ACDM background though). 

Aside from the LTB model, also a Szekeres model 
with no symmetries has been studied. This non- 
symmetric model is constructed by using the dipole 
functions defined hy S = 1, P = const, and 

choice leads to an overdensity peaking with a density 
contrast of approximately 4 at present time, 

t = to- 

2D present time density profiles of the LTB and 
Szekeres models are shown in figures 1 and 2. 


Using the above metric it is straight forward to 
obtain the geodesic equations in the spherical coordinate 
system: 



(23) 


Rie + R¥ + -k -k 0fc® -k 0/fc® - ^[P,^(fc’')2-k 

F,r(.k^f -k Pr-(P)^ + 2Q^rk^¥ + 2$,^P/U] = 0 

(24) 


Fk^ + Fk^ + Q¥ -k 0F - \[RAkl‘^+ 
PAk’^f + 2QA^¥ -k 2$,ePF] = 0 


( 25 ) 
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Pk'l’ + Pk^ + - ^[R^ik''f + 


(26) 


As with the Newtonian gauge metric, these equa¬ 
tions must be differentiated in order to obtain ODEs for 
, Pj ., pQ and fc'^ which are needed in the expression 
for Da- The resulting equations are: 

RtPPc + + Ftk^k%+ 

IPMk^f + P,tk‘^k% + + ^,tikiP+ 

k*KP + + Q,t{k<^,aP + Ppp] = 0 

R,c.P + Ri^riPP + P^ki) + RP^+ 

aA ’ ’ 

P{R,a0k^ + R,pPp + + ^,fik%)+ 

Hj^{k^P + k^pkP) + ^,o.kU 

+ Q{^K)+ ^ 28 ) 
k^k^P + e,aP + ek% - [^R,PP)P 

R.rPPc + IP^rPky + F,rPk% + + 

P,rPP^ + Q.rc^PP + eAk%P + PPA + 
^,rc.PP + ^Ak'^A" + PkA)] = 0 


Inserting the appropriate Christoffel symbols (given in 
appendix B) into equation (6), the following differential 
equation for the angular diameter distance is obtained: 


4^h^ = 2(fc* + P, + k% + kA+ 


dX 


,r 

Pr 


+ E^P + ?AP + ^P + ^P+ 
^ {R^p F R,rP + R,eP + R,A'^)+ 


R _ It! _ 
n p F 


0 


,0 


F R - ^ 20,t)+ 

^ p p 

P{^^-e - ^ 0,0 + ^F,r - 20 ,.) + k\-2QA+ 

A d> 1 $ 

fc^-P,-20,,)) + - ^_^_ 


F 


P 


(31) 




F 

2$ 0 + fc’’(|0.0 - + %P,r - 2$,.) 


+fc«(|pe - 2$,e) + fc^-2<i>,^ - |pe)) 

This equation is solved simultaneously with the 24 ODEs 
for the fc“’s and their derivatives. 


The last set of equations needed are again the null- 
condition and its partial derivatives. These equations 
are used when setting the initial conditions and for 
checking the accuracy of the code. 

In spherical coordinates, the null-condition is given 
by: 

fc“fc„ = k^kAo^p = -P{Pf + R{Pf+ 

F{Pf F P{kA‘^ + 2(dPP F 2^Pk‘^ = 0 

Taking the partial derivative of this equation one obtains: 


FpA F Fi^ik%) F k%k%) F k\F^^pA F F^pAA+ 
Fk% F P{Q,pA^ F Q^k^A + 0fc)a + 0.afc"+ 

0(^(fc:J + Ppk^A - [\R 9 APf + RePPo.+ 

\PfiAk*f + Pfik'^k*^ F e,0A^P F QAk%P+ 

PPA + ^,eA^P + ^AktP + k^PA] = 0 

(29) 


-2APk% F R,APf F 2RPPp + Fp{Pf+ 
2FpPp F Pp{Af F 2Pk*k% F 2epPpF 

2Q{PA^ F PAA + 2^pPA F 2^{PA^ F Pk^) = 0 

(33) 

The equations given above are shown in expanded 
form in appendix B together with comments on initial 
conditions. 


P-<.A F Pi^ik-^A + k'!;sk%) 
kAP-pA^ + P.pki) + Pki F P{^,apk^ + <i>,pki) 

U:^F^,o.PFA^{kA)FPAi) 
[IrMP)^ + R^PPp + QfA^’P + ©.^(fc'c.fc’' + k^P,A 

F ^AktP + ApA] = 
(30) 


C. Relation between the Szekeres and Newtonian 
gauge spacetimes 

In this section, the subscripts ”s 2 ” denotes Szekeres 
coordinates including LTB coordinates while the sub¬ 
script ”Zt6” is used for specifying LTB coordinates. The 
subscript ”ng" denotes Newtonian gauge coordinates 
and the subscript ’’//rw” is used for denoting coordi¬ 
nates on the unperturbed FLRW background. A tilde 






7 


will be used to denote coordinates of fiducial spacetime 
points. 

The numerical values of coordinates are of no physical 
value. Instead, the physically relevant quantities are 
proper distances, and these are determined by the 
metric. In order for two spacetimes to be considered 
equivalent, a map between their coordinates must 
thus be constructed such that proper distances are 
equal in the two spacetimes. It could seem, that the 
appropriate point identification map in this case is 
between the Szekeres and Newtonian gauge spacetimes 
and this would indeed also lead to an interesting study. 
However, the recipe of [17] is for going between the 
Newtonian gauge metric and an N-body simulation with 
its underlying FLRW metric. The intent here is to see 
how well the recipe describes the ’’true” underlying 
relativistic spacetimes corresponding to density and 
velocity fields obtained from N-body simulations. This 
is done by studying the Newtonian gauge metric’s 
ability to accurately describe light propagation. Such a 
study requires the knowledge of the ’’true” underlying 
spacetime of the considered N-body fields. Thus, mock 
N-body data is constructed by mapping Szekeres models 
into their background FLRW models. The mapped 
Szekeres model’s density profile and the velocity field 
obtained through the mapping comprise mock N-body 
data. As in regular perturbation theory, the Newtonian 
gauge metric is then assigned the same coordinate 
system as that of the FLRW background of the ” N-body 
data”. 

Both the velocity and density fields of LTB models 
can be reproduced by N-body simulations when using 
initial conditions based on the maps described below. 
As mentioned in the introductory section, it is still 
work in progress to show that this is also the case for 
non-symmetric Szekeres models. Following the standard 
consensus, it is here assumed that Newtonian N-body 
simulations accurately reproduce structure formation 
in accordance with general relativity. In particular, 
it is thus assumed that non-symmetric quasi-spherical 
Szekeres models are reproducible by Newtonian N-body 
simulations. 

Following [20], the point identification map be¬ 
tween FLRW and LTB coordinates is defined by the 
requirement of equal proper radial distances in the 
two spacetimes. Letting denote components of 
the metric tensor, the identification between comoving 
coordinates in the two spacetimes is thus given by the 


four equations: 

t := tlth — tflrw 

^ ^Itb — ^flrw 

4^ •— 4^ltb — 4^flrw 
j-fitb rfjirui 

dPrptb I — I drflrwy/9rr,flrw 

Jo Jo 

(34) 

These equations yield a mapping (t,ffirw,0,(j)) —>■ 

{i,rith,0,(})). Note that since the big bang time in 
equation (18) is zero, the time coordinates of the two 
spacetimes are the same. Note also, that the spherical 
symmetry about the origin of the LTB model implies 
that the angular coordinates of the LTB and FLRW 
metric are the same. 

The non-symmetric Szekeres metric is non-diagonal 
in spherical coordinates while the FLRW metric is 
diagonal. A point identification map in spherical coor¬ 
dinates does thus not capture the complete anisotropy 
of the Szekeres model. As shown in equation (15) the 
general Szekeres metric is diagonal in stereographic 
coordinates defined by equation (14). The appropriate 
point identification map between the non-symmetric 
Szekeres spacetime and the FLRW spacetime is thus 
defined in stereographic coordinates: 

i .— isz — iflrw 
P ■— Psz ~ Pflrw 
Q Qsz — Pflrw 

prsz nTflr-w 

PPr^sz ■ / szy/9rr,sz / flrw yjQrr^flrw 

JO Jo 

(35) 

In this equation, grr,sz is the rr-component of the Szek¬ 
eres metric in stereographic coordinates while grr,itb was 
the rr-component of the LTB metric in spherical coor¬ 
dinates. The rr-component of the FLRW metric is the 
same in spherical and stereographic coordinates as the 
flat FLRW metric in stereographic coordinates is given 
by: 

ds^ = —C^dt'fir^+a'^ + ^2 ^^p}lrw + dqjirw)^ 

(36) 

The function E in this line element is defined hy E = 
5 {pflrw + d}irw + l) which in spherical coordinates cor¬ 
responds to E = 2 sin^( 6 i/i — Ji) • 3"be stereographic FLRW 
coordinates are related to the spherical FLRW coordi¬ 
nates by the transformation: 


Pflrw — cot (^0firw/^') COs(0J/7-^^) 
Qflrw ~ cot {Ofirw/^) ^^^iPflrw) 


( 37 ) 








Combining the stereographic point identification 
map with the coordinate transformations of equations 
(14) and (37), a point identification map between the 
spherical coordinate systems of the Szekeres and FLRW 
spacetimes is obtained. The stereographic coordinates of 
the Szekeres model are related to the angular coordinates 
through an r-dependence but this is not the case for the 
FLRW metric. This difference implies that the angular 
coordinates in the two spacetimes will not be identical. 

Other point identification maps than the ones given 
above have been studied. In particular, a map requiring 
equal proper distances along the (p, g)-coordinate direc¬ 
tions instead of Psz = Pfirw,<lsz = Qfirw has been studied 
along with two maps in spherical coordinates. For the 
non-symmetric Szekeres model, the maps do not yield 
the same results indicating a lack of exact compatibility 
between the FLRW and general Szekeres spacetimes. 
The map shown in equation (35) is used because it yields 
the best reproduction of the non-symmetric Szekeres 
spacetime. 


a. Peculiar velocities and In the Szekeres model, 
the observer is comoving and the redshift is given by 
I + z = ^7^, with a subscript 0 referring to evalu- 


ation at the position of the present time observer. In 
the Newtonian gauge spacetime, dust is generally not co¬ 
moving. Following [17], the appropriate peculiar velocity 
field is proportional to (I,u*) with u® = Vi the local co¬ 
moving spatial dust motion on the EdS background/in 
the N-body simulation. The velocity field is normal¬ 
ized such that u^Uu = —(? so the appropriate field is 




The 


velocity field is needed in order to obtain the redshift 
along the geodesics using the general formula 1 + z = 
(fc°M°)o ’ ^ subscript e indicating the spacetime po¬ 

sition of emission i.e. = fc“(A) etc. 

The comoving peculiar radial velocity Vr is computed 
by taking the time derivative of the expression for the 
proper distance dpr'. 


= ^{arfirw) = a^trfirw + avr (38) 

The proper distance on the left hand side is computed 
in the Szekeres/LTB spacetime at the appropriately 
mapped spacetime point. 

The angular velocities vanish in the LTB model 
because of the spherical symmetry about the origin. The 
angular velocities of the non-symmetric Szekeres model 
do not vanish identically, but they are small and are 
only used in the formula for the redshift in which they 
are suppressed by k^g^kf^g- The angular velocities can 
thus be neglected. 


The potential ipng is needed in order to solve the 


geodesic equations and for obtaining Da in the New¬ 
tonian gauge. The potential is obtained by the usual 
Poisson equation W^ipng = dpng- This equation 

is solved in Fourier space (using fftw3^) on a grid and 
the value of tpng at a specific point is then obtained 
using quadri-linear interpolation. This method should 
approximately mimic how one would work with the 
potential from an N-body simulation. 

The overdensity is defined as 5p := p Szekeres — PEdS 
and is obtained at spacetime points in the Szekeres 
metric. The overdensity needed for obtaining ipng is 
the corresponding overdensity of the mock N-body 
simulation so a mapping into the Newtonian gauge 
spacetime is necessary. The overdensity 6png at a New- 
tonian gauge spacetime point {ing,rng,Ong,}ng) = 
flrw) ^ firm : ^ flrw) ^ flrw') tllUS COmputed. aS 

^Pngij^ngi^ngt ^ng-) ^ng^ — ^PSzekeres itsz: ksz: dsz: ^sz)j 

where {tsz,fsz,dsz,4’sz) is the Szekeres spacetime point 
corresponding to {tngifng,dng,(png) according to the 
point identification map defined by equation (34) or (35). 
The corresponding potential in the Newtonian gauge 
spacetime, ipng, is then computed from the Poisson 
equation in the Newtonian gauge spacetime. 


D. Gauge transformation of the angular diameter 
distance 

In the previous subsection, it was discussed how 
to construct a point identification map between the 
Szekeres model and its FLRW background in order to 
construct mock N-body data. This map is used for 
obtaining the metric potential ip and peculiar velocities 
of the Newtonian gauge metric so that this metric de¬ 
scribes the Szekeres spacetime. Another issue regarding 
the use of the two spacetimes is that they do not corre¬ 
spond to the same spacetime slicing since the Szekeres 
metric is given in the comoving synchronous slicing of 
spacetime. In order to compare Da{z) obtained from 
the Newtonian gauge equations with the exact result, it 
must thus in principle undergo a gauge transformation 
to the comoving synchronous gauge. 

It is well known from linear relativistic perturbation 
theory that the differences in observable quantities in 
different gauges are negligible well within the horizon. 
The structures studied here have a present time diameter 
of approximately IGOMpc so the gauge transformation 
should not be necessary. A formal gauge transformation 
is performed anyway as a precaution. 

Redshift perturbations are gauge dependent, but 
the redshift itself is a scalar i.e. coordinate invariant. 
Since the value of the exact redshift is computed here. 


® http://www.fftw.org/ 
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and not just the perturbation, no transformation of 
the redshift is thus necessary. The angular diameter 
distance is not gauge invariant; it is the fraction of two 
infinitesimal areas and areas change under coordinate 
transformations. 

As will become apparent below, only the transforma¬ 
tion of the time coordinate is needed to obtain a gauge 
transformation of the angular diameter distance. The 
transformation law between the time coordinates in the 
two coordinate systems is found by using the general 
coordinate transformation equation = 9ai3§^§^- 
Inserting —as the time-time component of the Syn¬ 
chronous gauge metric tensor into the left hand side 
of this equation and inserting the Newtonian gauge 
metric into the right hand side, yields the approximate 
transformation law: 

ntN 

is / i^dt (39) 

Jo 

In this equation, tjv and ts are the time coordinates in 
the Newtonian and Synchronous gauges respectively. 

The approximate result is the same as what one would 
obtain by using regular gauge transformation laws for 
going between the Newtonian gauge and a synchronous 
gauge (see e.g. section 5.3 of [30]). By using the first 
order Euler equation in the spherically symmetric limit 
it is seen that = —a f^Vr(t]^,r')dr' which 

shows that the approximate result obtained here is also 
in agreement with that found in [31]. 

Using the simplihed notation Xg = x'^ + for 

the gauge transformation of x°‘ = a Taylor 

expansion reveals the gauge transformation of the 
angular diameter distance: 

+ (5A«,0).ario 

An over-bar is used for denoting background quantities 
and coordinates without a subscripted N or S are 
background coordinates. As before, N and S denote 
coordinates in the Newtonian and comoving synchronous 
gauge respectively, and e denotes point of emission while 
0 denotes point of observation. 


Using 


2cH, 


-1 


the 

a(*e) 

a(to) 


law becomes: 


Mattig 

a(t,)3/2' 
a(to)3/U ’ 


relation [32], Da = 
the gauge transformation 


Da,s{ 


e,S’ ^0,S 


j) ~ Da,n{ 


e,N ) ^0,N 


r) + 


2cH-^ 


(u) 

a{to) 


3 yja{te)a^t{te) 
2 aito^D 


e\e 


ce‘lo 


(41) 


Angular diameter distance along radial LTB ray 



FIG. 3. Angular diameter distance along a radial geodesic in 
the LTB model with a central observer. The line with the 
legend ”LTB ray” is the angnlar diameter distance along a 
geodesic computed with the ODE’s of the exact LTB met¬ 
ric. The two other lines are the angular diameter distances 
along the "equivalent” geodesic according to the Newtonian 
gauge metric in the Newtonian gauge and in the synchronous 
gauge respectively - these two lines are indistinguishable as 
expected. They are also indistinguishable from the first men¬ 
tioned line indicating that the Newtonian gauge metric ade¬ 
quately describes the LTB geodesic. 


non-symmetric Szekeres models described by the New¬ 
tonian gauge metric. Afterwards, the exact Da{z) re¬ 
lation is obtained by solving the set of ODEs given in 
section IIIB. In order to compare the geodesics and the 
corresponding distance-redshift relations obtained with 
the two different sets of ODEs, the geodesics must of 
course be initialized equivalently. Since the equations are 
solved backwards in time, this implies that the observer 
position and line of sight must be mapped from the New¬ 
tonian gauge coordinate system to Szekeres coordinates. 
The position of the observer is mapped using the maps 
of equation (34) or (35). The line of sight of the observer 
is mapped by mapping the Newtonian gauge geodesic 
into Szekeres spacetime and computing a two point finite 
difference along the beginning of this geodesic. The two- 
point finite difference constitutes the initial conditions 
of (fc’", fc®, fc"^) and the initial condition of fc* is obtained 
through the null-condition. A geodesic which is initial¬ 
ized as radial at the position of a central observer is radial 
in both spacetimes. No mapping of fc'’, fc®, is necessary 
in this case since the results are frequency independent 
and fc* and D' uniquely determine each other through the 
null-condition. The initial direction of the ray is deter¬ 
mined by a set of angular coordinates {9, (jj) which must 
be mapped though. 


A. Geodesics in the LTB spacetime 


IV. RESULTS 

The system of ODEs presented in section II are used 
to obtain Da{z) along single geodesics in the LTB and 


In figures 3 and 5, Da{z) is shown for a radial and a 
non-radial geodesic in the LTB model respectively. The 
radial geodesic corresponds to a central observer, while 
the observer is placed outside the void in the non-radial 
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Density along radial LTB ray Density along nonradial LTB ray 




FIG. 4. Density along a radial geodesic in the LTB model with 
a central observer. The density is plotted both according to 
the Newtonian gauge geodesic and the Szekeres geodesic, but 
the two densities are the same implying that the two geodesics 
pass through equivalent portions of spacetime. 


Angular diameter distance along nonradial LTB ray 



FIG. 5. Angular diameter distance along a non-radial 
geodesic in the LTB model with an observer placed outside 
the void at r = 1.2 in Newtonian gauge coordinates (comoving 
coordinates are normalized at t = tis in the units of O.lMpc). 
The geodesic is initialized in the spatial direction determined 
by (fco, ^ 0 i ^ 0 ) = (—0.01, 0.001,0) in the LTB spacetime. The 
line with the legend ”LTB ray” is the angular diameter dis¬ 
tance along a geodesic computed with the ODE’s of the exact 
LTB metric. The two other lines are the angular diameter 
distances along the ’’equivalent” geodesic according to the 
Newtonian gauge metric in the Newtonian gauge and in the 
synchronous gauge - as expected, these two lines are indis¬ 
tinguishable. They are also indistinguishable from the first 
mentioned line indicating that the Newtonian gauge metric 
adequately describes this non-radial LTB geodesic. 


case. For the radial ray, there are only negligible 
differences between the curves obtained from solving the 
ODE system based on the Newtonian gauge metric and 
that based on the exact LTB metric. For the non-radial 
ray, there is a slightly more noticeable difference between 
the distance-redshift relation obtained with the exact 
metric and that obtained with the Newtonian gauge 
metric. This difference is presumably due to precision 
errors occurring when computing the initial values of 
,k^, k^ in one spacetime from their values in the other 
spacetime. 


FIG. 6. Density along a non-radial geodesic in the LTB model 
with the observer placed at comoving r = 1.2 in Newtonian 
gauge coordinates (comoving coordinates are normalized at 
t = tia in the units of O.lMpc). The observer is looking in 
the direction determined by (fco,feQ,fcQ) = (—0.01,0.001,0) 
in the LTB spacetime. The density is shown along both the 
Szekeres and Newtonian gauge geodesics. The two densities 
are the same though, indicating that the Newtonian gauge 
adequately reproduces the LTB geodesic. 

In the figures, Da{z) is also shown after a transforma¬ 
tion to the comoving synchronous gauge. Clearly, this 
transformation is completely obsolete as was expected. 

In figures 4 and 6 the density profiles along the 
rays are shown ^ . These are the same along the exact 
LTB rays and the Newtonian gauge rays which shows 
that the rays in the two spacetimes follow equivalent 
spacetime paths. 

It is not surprising that the distance-redshift relation 
of the LTB model is reproduced by the perturbed FLRW 
model in the Newtonian gauge, as it has earlier been 
shown that the Newtonian gauge describes the LTB 
spacetimes well [31, 33, 34]. 

B. Geodesics in non-symmetric Szekeres 
spacetimes 

In figure 7, the distance-redshift relation along a 
geodesic in the non-symmetric Szekeres model is shown. 
The angular diameter distance is not as precisely 


^ Note that in the abridged dictionary of [17], the density used 
for computing the potential is not the same as the density of 
the Newtonian gauge spacetime. Using the notation of [17], 
<5 = <5jv — 4 .n-Qp^ ds°-‘^ + O',tip,t), where 5 is the density 

contrast corresponding to the Newtonian gauge spacetime while 
5jv is the density contrast of the N-body simulation and i/ijv the 
corresponding potential. V’iv is the potential perturbation ap¬ 
pearing in the metric, while S describes the density field of the 
spacetime. In the notation used here, corresponds to Sng 
which was used to compute the potential. The difference be¬ 
tween the two density fields of the dictionary is insignificant for 
the models studied here. 
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Angular diameter distance along Szekeres ray 



Polar angle along ray in non-symmetric Szekeres model 



FIG. 7. Angular diameter distance along a geodesic in 
the non-symmetric Szekeres model with a central observer. 
The ray is initialized as radial in the direction specified by 
(0o,<^o) = (7r/3,7r/3) in Newtonian gauge spacetime coordi¬ 
nates. The line with the legend ’’Szekeres ray” is the angular 
diameter distance along a geodesic computed with the ODE’s 
of the exact Szekeres metric. The two other lines are the 
angular diameter distances along the ’’equivalent” geodesic 
according to the Newtonian gauge metric in the Newtonian 
gauge and in the synchronous gauge. A close-up of the area 
where the exact and Newtonian gauge results begin to be dis¬ 
tinguishable is included. 


Density along Szekeres ray 



FIG. 9. Polar angle along the geodesic in non-symmetric Szek¬ 
eres model. The angular changes are shown along both the 
exact ray according to the Szekeres metric (’’Szekeres ray”) 
and along the ray according to the Newtonian gauge metric 
(’’Newtonian gauge ray”). For the Newtonian gauge geodesic, 
the angular changes are shown in both Szekeres and Newto¬ 
nian gauge coordinates, denoted by ”sz” and ”ng” respec¬ 
tively. 


consistent with its moving through a larger overdensity 
at the void edge and implies that the Newtonian gauge 
ray and the exact ray do not move along exactly equiva¬ 
lent spacetime paths. Similar results are obtained when 
studying rays in Szekeres models with less anisotropy 
and smaller density contrasts, with the disagreement 
between the Newtonian gauge and exact Szekeres results 
becoming less as the level of anisotropy and overdensity 
are decreased. Models with more anisotropy have not 
been studied. 

As with the LTb rays, the gauge dependence of 
the distance-redshift relation along geodesics in the non- 
symmetric Szekeres models is completely insignificant. 


FIG. 8. Density along a geodesic in the non-symmetric Szek¬ 
eres model with a central observer. The ray is initialized at 
a central observer as a radial ray in the direction specified by 
{Oo,<j>Q) = (7r/3,7r/3) in Newtonian gauge spacetime coordi¬ 
nates. 


reproduced by the Newtonian gauge metric as in the 
LTB case, but the reproduction is still very good. The 
difference between the Da{z) relations along the two 
rays is consistent with the Newtonian gauge ray moving 
through slightly larger overdensity than the exact ray. 
This is in fact the case, as can be seen in figure 8 which 
shows the overdensities along the two geodesics. In 
figure 9 the polar angle along the two geodesics is shown. 
The geodesic is not bent in the Newtonian gauge metric, 
but when mapping the angular coordinates into the 
Szekeres coordinates, the geodesic is seen to correspond 
to a non-radial Szekeres ray. The ray is not bent exactly 
as the exact ray though. 

The slight under-bending of the ray inside the void is 


V. CONCLUSIONS 

The Newtonian gauge metric corresponding to the 
recipe of [17] was used to derive a set of ODEs that can 
be solved to obtain distance-redshift relations in the 
corresponding spacetime. The equations were used to 
obtain the angular diameter distance as a function of 
the redshift in mock N-body simulations based on quasi- 
spherical Szekeres models. The validity of the recipe 
was then estimated by comparing with redshift-distance 
relations obtained by using the exact Szekeres metric. 
In the spherically symmetric case, the distance-redshift 
relation obtained with the Newtonian gauge metric and 
the exact metric of the model agreed. This is not sur¬ 
prising since others have earlier shown that LTB models 
after a gauge transformation at least approximately can 
be described as spherical perturbations on an FLRW 
background. 

The results obtained here emphasize that it is futile 
to try to refute results obtained with LTB models by 
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claiming that they are unreliable because of spurious 
gauge modes in the synchronous gauge; it is a well known 
result from linear relativistic perturbation theory that 
observables have negligible gauge dependence on scales 
below the horizon. Here it has been shown explicitly 
that this is specifically the case for distance-redshift 
relations in quasi-spherical Szekeres models - even when 
the density contrast moves into the non-linear regime. 

For the non-radial geodesic in the non-symmetric 
Szekeres model, there are small but definite differences 
between the geodesic paths and distance-redshift rela¬ 
tions obtained with the two metrics. The ray appears 
not to being bent enough inside the void which leads it 
to move through a slightly larger density contrast than 
the exact ray. This again induces a slight difference in 
the distance-redshift relations obtained along the exact 
and Newtonian gauge geodesics. 

Versions of the studied recipe are the typical basis 
for ray tracing schemes. Standard ray tracing methods 
do not include treating the Newtonian gauge metric 
as exact though, and typically involve several different 
approximations. For example, a ray is typically not 
traced by the actual Newtonian gauge geodesics, but is 
traced in the background FLRW metric, possibly being 
bent at a finite number of lens planes. This may lead to 
results that are in less agreement with the exact results 
than what was obtained here. In addition, the small 
discrepancies between the exact and Newtonian gauge 
rays found here may be cumulative though it seems just 
as likely that the effects will cancel each other if the rays 
move through several structures instead of just the one 
they moved through in this work. 

To study how accurate actual ray tracing techniques 
are in reproducing shear, magnification etc. in exact 
spacetimes, a study of quasi-spherical swiss cheese 
Szekeres models is under way. If any shortcomings 
appear and especially if they are statistically robust, 
more accurate ray tracing techniques may need to be de¬ 
veloped in order to meet the increasingly high precision 
of observations. In such a case, exact inhomogeneous 
and anisotropic cosmologically relevant models seem like 
a good tool for general developments and tests of more 
accurate ray tracing methods e.g by studying the use of 
more advanced metric approximations such as one based 
on a post-Friedmann expansion (see e.g. [14]) or the 
’’Oxford” dictionary in [17]. 
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Appendix A: DDEs for the perturbed FLRW metric 
in the Newtonian gauge (Spherical coordinates) 

In this appendix the set of ODEs used to obtain 
Da{z) for perturbed FLRW metrics in the Newtonian 
gauge in spherical coordinates are given. 

The four geodesic equations are: 

c^(l -I- 2ip)k* = —2c^ipk* + 

(1 — 2i/')aa_t][(fc’')^-I--I-sin^(0)(/c‘^)^] 

(1 - 2il))a^k'' = 2a^'^¥- 

2(1 — 2'>p)aa,tk*¥ — — a^i/',r.(fc’’)^-|- (A2) 

a?r[(l — 2'ip) — + sin^(0)(fc'^)^] 

(1 - 2i))a^r^k^ = 2'iPa‘^r‘^k^ - 2(1 - 2il})ark\ra,tk'^ + a¥)- 

+ r2(fc«)2 + An\e){k*f]+ 

(1 — 2il))a^¥ sin(0) cos(6>)(fc‘^)^ 
(A3) 

(1 — 2'!/)a^r^ sh¥{9)k'^ = 2'ipa‘^r'^ sh¥{9)k'^— 

2(1 — 2'il;)k’^arsm{9)[a^tk^rsm{9) + asm{9)¥ + ar cos{9)k^] — 

cV 0(fc‘)' - ^ + r\k^f + sin\9){k^f] 

(A4) 

Taking the partial temporal derivative of equation (Al) 
one obtains: 

c^(l -I- 2i/)) —= —c^(l -I- 2'^l})¥^k^^— 

— 2(?¥{'>p^tpk^ + '^’.pk^t) ~ 

+ 2c'^'tfj^tk*k*t + 2 [V',ta^ ~ (1 ~ 2ij;)aa^t][¥¥t~k 

r'^k^k^^ + s\t?{ 9)k‘^k'^f\ + -I- 4'if:^taa,t— 

(1 - 2ip){a\ -I- aa^tt)][{¥)^ + r‘^{k^f' + sin^(6>)(fc‘^)^] 

(A5) 

The partial r, 9 and (j) derivatives of Al corresponds to 
the following three equations: 

c^(l -I- 2ij;)-^¥ = —c^(l -I- 2ip)¥gk^ — 2c^-!/’rfc* — 
dX ' ’ 

2c^¥[(l)^rpk^ + ''P,pk^r\ ~ ‘2c^i’k*r + C^'4’,tr{k*)‘^ + 
2c^tp^t¥^k* + 2['!/_ta^ — (1 — 24>)aaA][¥ ¥r + r{k^)‘^+ 
k^k^., + r sin'^{9){k’^f+r^ sin^{9)k’^k'l’^] + 

2^ raa i][(F)" + r^{ky + sm^{9){k’^f] 

(A6) 
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c^(l + 2^)-^k*g = —c^(l + 2'ip)kKk^g — 2c^i/;^efc* — 
2c^/c‘[V',aefc“ + V',afc^] - 2c^ijjk*g + c^i’^t0(k^f + 

2c^'ip^tk*k'^g + 2['(/',ta^ ~ (1 ~ 2ip)aa^t][k''k^g + r‘^k^k^g + 
sm{6) cos{d){k'^)‘^ + sin^(0)fc'^fc g] + [^/’,tea^ + 

2^,gaaMkl^ + r^k^ + sin^{e){k^f] 

(A7) 

c^(l + 2-(/')^fc*0 = -c^(l + 2'ilj)k*pk^^ - 2c^V'.(/>^*- 

2c^fc*[^/’,a0fc“ + - 2c^i’k*4, + c^'^l^,t(i,{k*Y + 

2c^^,tfc‘fc*^ + 2[-!/',ta^ - (1 - 2V')aa,t][fc’'F^+ 
r'^k^k^^ + sin^(0)A:'^A:|^] + [-!/',t0a^+ 

2^ ^aa t]Wf + r^iky + sin\e){k'^f] 

(AS) 

The partial derivatives of equation (A2) are: 

(1 — 2tp)a^ ^k'^^ = —[2(1 — 2ilj)aa^t — 2V’,ia^]fc’'— 
uA 

(1 — 2 'i/))a^fc(gA:^ + 2[2aa^t4’ + o,'^i'4’,i3k^t+ 

V’,t/3fc^)]fc’’ + ‘^a^'>pk''t - 2[(1 - 2'!/')(a,tta+ 
a^t) — 2ijj^taa,t]k*^k'^ — 2(1 — 2tjj)aa^t{k*tk''+ 
fc*A;[)) — (?il)^tr{k^Y ~ 2.c^'ip,rk*k*t — [2aa^t'<P,r+ 
o^ip^tAik'")^ — 2a^V',r-A’'fc[t + 2 aa_tr[(l — 2 ij}) — rip ^j.][{k^y+ 
sin^(0)(A:'^)^] + 2a^r[(l — 2ip) — rtjj ,r][k^ k^t+ 

sin^(0)fc'^fc‘^] + a'^r[—2'ip t — rip tr][(fc®)^ + sin^(0)(fc''^)^] 

(A9) 

(1 - 2iP)a^-^¥^ = -(1 - 2iP)a^¥pkP^ + 2iP^ra^¥+ 
clA ’ ’ ’ 

2a^ipkp^ + 2a^¥[ip^rpk^ + ^,/3fc[^] + ^4’, raa^tk^¥— 

2(1 — 2ip)aa^t[k*rk’" + k*¥^] — c^ip^rr{k*Y ~ 2,c^ip,rk*k*j.— 

a^ip^rr{¥Y ~ 2a^ip^r¥kp^. + a^[(l — 2ip) — rip^r\[{k^Y+ 

sin^(0)(A:‘^)^] + a^r[—Si/'.r — »'V',rr] [(fc^)^+ 

sin^(0)(fc‘^)^] + 2ra^[(l — 2ip) — rip^r\[k^k^r + siv? {9)k'^ k%] 

(AlO) 

(1 — 2ip)a'^ -^¥g = —(1 — 2ip)a^¥gk^g+ 

2a^ip^g¥ + 2a'^¥[ip^i3gk^ + ip,pk^g] + 2a^ipkPg — 
2(1 — 2ip)aa^t[k*^ek^ + k*¥g] + Aipgaa^tk*¥— 
<?4',re{k*y — 2¥ip^rk*k*g — a^ip^re{¥Y ~ ‘2a‘^ip^r¥kpg+ 
a‘^r[-2ip^g - rip^r0][{k^)'^ + sin^(0)(fc‘^)^] + 2a^r[(l - 2ip)- 

rip r][k^k% + sin(0) cos(0)(fc'^)^ + sin^(0)fc''^fc'^] 

(All) 


(1 - 2iP)a^ — kp^ = -(1 - 2iP)a^kPpk^^+ 
2a^ip^^¥ + 2a^¥[ip^p^k^ + ip,pk^^^ + 2a^ipkp^- 
2(1 — 2ip)aa^t{¥^¥ + k*kp^) + Aip ^tk* ¥ — 

c^ip^r<p{k*)'^ — 2¥ip^r-k*¥^ — a^ip^r<t>{k''Y ~ 2a^ip^rk^¥^+ 

a^r[- 2 V ',0 - rip^rc/>][ik'^)^ + sin^(6>)(fc‘^)^] + 

2a^r[(l — 2ip) — rip^r][k^kP^ + sin^(0)fc'^fc‘^] 

(A12) 


The partial derivatives of equation (A3) are: 


(1 - 2ip)a^r'^-^k'pt = -(1 - 2ip)a^r‘^k‘^igk^t- 

¥k^[—2ip^tO^ + 2(1 — 2ip)aa^t] + 2r^A:®[2aa^t'0+ 
o?{ip^ptk^ + V’,/3^Jt)] + ‘2a^r'^ipk^^ — 2[(1 — 2ip)a^trk^ — 
2arip^tk^ + (1 — 2ipi)ark^P\\ra ,tk*' + a¥] — 
2(1 — 2ip)ark^[ra^ttk* + ra^tk*t + ri,t¥ + a¥^] — 
c^'fp,te{k*)^ - 2¥ip^g¥kP - [2aa^t'<P^g + a^^/>_(e][(F)^+ 
r2(fc®)2 + sin2(0)(A:<^)2] - 2^iP^g[¥¥^ + r‘^k*^k^^^+ 

¥ s\r?{6)k'^k^P\ + 2r^ sin(0) cos(6*)[(l — 2ip)aa^t{k'^)^ — 

ip,t¥{k'^f + {l-2iP)¥k'^k^t] 

(A13) 


(1 - 2ip)¥¥^k% = -(1 - 2ip)¥¥k%k^,- 
a^fc®[2r(l — 2ip) — 2¥ip¥ + 2¥k^[2rip+ 
¥{ip^rpk^ + i’,pk^r)] + 2.¥r^ipkPr- 
2a[k^{l — 2ip) — 2rip^rk^ + (1 ~ 2ip)rkP^][a^trk* + a¥] — 
2(1 — 2ip)ark^[a^tk* + ra^tk*^. + a¥^] — ¥ip^gr{k*)‘^ — 

2¥ip^g¥¥^ - ¥ip^gr[{¥f + ¥{k<^f + ¥ si¥{e){k*f]- 

2¥ipfi[¥¥^ + r{k^)'^ + ¥k^kPj.+ 
rsi¥{9){k‘^)‘^ + ¥ si¥(9)k'^k'p.] + 2¥ cos(9) sin(9) 

[r(l - 2ip){k'^f - ¥ip¥k‘^f + (1 - ‘2iP)¥k^ki] 

(A14) 
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(1 — 2ip)a^r‘^-^k^g = —(1 — 2'ip)a^r'^k%k^g+ 

+ 2a^r'^k^[il)^pgk^ +'0,^fcg] + 
2'^a^r'^k^g — 2or[(l — 2ip)k^g — 2ilj^gk^][ra^tk* + ak^] — 
2ar{l — 2'tfj)k^[ra^tk*g + ak^g] — c^'ijj^gg{k*)‘^ — 
20^11)fik^k^g - a^il}fig[{k''Y + r^(fc®)^+ 
sin^(0)(fc‘^)^] — 2a?'ij)^g[k''k'^g + r'^k^k^g+ 
sin(0) cos(0)(fc'^)^ + sin^ iO)k'f>k^] + 
a‘^r^{k'^)^[{l — 2^)(cos^(0) — sin^(0)) — 2'0,e sin(0) cos(0)] + 

2(1 — 2i/))a^r^ sm{6) cos{9)k'^k^g 
(A15) ’ 


(1 - 2^)a^r‘^^k'^^ = -(1 - 2i})a^r'^k^pk^^+ 

+ 2a^r^^/>fc®0 + 2a^ k^ [li) ^p^k^ + i’^pk^^] — 
2[(1 — 2%lj)ark^^ — 2ar'ip^^k^][ra^tk^ + ak''] — 
2ar{l - 2V')fc®[ra,tfc‘^ + ak^^] - c^il)^g^(k*Y- 
2cV.e/c‘fc*^ - aV.e0[(fc")2 + r2(fc®)2 + r2 sin\e){k^f]- 
2a^ipfi[k"'k'^^ + sin^(0)fc'^fc^]- 

sin(0) cos(0)(A:'^)^ + 2a^r^(l — 2^) sin(0) cos(0)fc'^fc‘^ 

(A16) 

The four partial derivatives of equation (A4) are: 

(1 — 2^/>)o^r^ sin^(0)^fc^ = 

CLA 

— (1 — 2^)a^r^ sin^(0)A:‘^/c|^ — 

2r^ sin^(0)afc‘'^[—-0,(0 + (1 — 2^/^)a,t] + 

2r^ sin^(0)[a^fc‘^('!/:^Q,(fc“+ 
tp^ak'^t) + ‘2'aa^tk^'ip + a^tpk'll] — 
2rsin(0)[(l — 2'ip)a^tk'^ + (1 — 2ip)ak‘l'f — 
2i/)_tafc‘^][a_t/c*rsin(0) + asin(0)fc’’+ 
arcos{9)k^] — 2(1 — 2i/')ar sin(0)fc'^[a_ttfc‘r sin(0)+ 

a.tfc^jT sin(d) + o_t sin(d)/c'' + asin(0)A:((+ [ 

a^tf cos{9)k^ + arcos(d)fc®] — (?'ip^pt{k*'Y — < 

2c^V'.0^‘fc*t - + 2aa,tV’,0][(^’')^+ 

r2(fc®)" + r" sin2(0)(fc'^)2] - 2a^pj^4k^kPt+ 

r^k^k\ + sin^(0)fc'^fc‘^] 

(’a17) 


sin^(0)(l — 2^/>)^fc(^ = 
ctA ’ 

—a^r^ sin^(0)(l — 2'0)fc'^fc(^— 

2a^r sin^(0)fc'^[l — 2'!/) — ■0,,?']+ 

2a^ sin^(0)fc'^[2r^/> + r^(V',r/3^^+ 
i/;,/ 3 fc(^)] + 2ipa'^r‘^ sin^(0)fc(^— 
2asin(0)[(l — 2'ip)k'^ — 2'iljrfk^+ 

(1 — 2^/))rA:|^][a_tfc‘r sin(0) + asin(d)A:’’+ 

ar cos{9)k^] — 2(1 — 2'ip)k‘^ar sin{9) (^1^) 
[a^tk^j.r sin(0) + a^tk* sin(0) + a sin(0)/c((,+ 
a cos{9)k^ + ar cos(0)fc® ] — 
c^'0,0r(fc‘)^ - 2c^tp^^k*k*^ - ip,r4,a^[ik^)‘^ + 

r^{ky + 

sin2(0)(A:<^)2]- 

2o?p)^^[k''kPj. + r{k^Y + r'^k^k^r+ 
rsin^(0)(fc‘^)^ + sin^(0)fc‘^/c|^] 


(1 — 2^/>)a^r^ sin^(0)^fc'^ = 

— (1 — 2ip)a^r'^ sm^(9)k'^i^k^g — 2a^r‘^ sm{9)k'^ 
[—tp^e sin(0) + (1 — 2'ip) cos(0)] + 2a^r^k’^ 
[2^sin(0) cos(0) + s\T?{9){'ip^pgk^ + 'ip,pk^g)]+ 
2ipa?r‘^ s\T?{9)kP^g — 2ar[(l — 2'ip) sm{9)k'^g + 
(1 — 2'ip)k'^ cos{9) — 2'ip^gk'^ sm{9)][a^tk*r sm{9) + 
asm{9)k'~ + ar cos{9)k^] — 2(1 — 2ip)k‘^ar s\n{9) 
[a_t/c‘g sin(d)r + a^j/c* cos(d)r + acos{9)k'^ + asm{9)kpg + 
ar cos{9)k^g — ar sm{9)k^] — c^'ip^^g{k*)'^ — 2c^tp^^k*k*g — 

ip,4,ea'^[{kJ')‘^ + r^(fc®)^ + 
sin^(0)(fc‘^)^] — 2a'^'ip^p[k^kpg + r‘^k^k^g + 

sin^(0)fc'^fc g + sin(0) cos(0)(fc‘^)^] 

(A19) 
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(1 - 2'tp)a'^r‘^ sin2(6»)^fc^^ = 

— (1 — 2'tp)a^r'^ sin^(0)/c^/c|^+ 
2a?r^ shi^{0)ip ,j,k'^ + 2a?r^ sin^(0)'!/>/c‘^+ 
2a^r^ sir?(B)k'*‘[(l)^p^k^ + 

2arsin(0)[(l — 2'(/i)/c‘^ — 2'ilj^^k'^] 
[a_tfc‘r sin(0) + asm{9)k'^ + ar cos{9)k^] — 

2(1 — 2i/))fc‘^arsin(0)[a_tr sin(0)A:*0+ 
asiii(6i)fcC^ + arcos(6»)fc®^] - c^'ii)^^^{k*f- 
2c^il)^^k*k^^^ - + r‘^{k^Y + 

sin^{e){k^f] - 2a^4,,^[k^k^^ + r\‘>k%+ 

sin^(0)fc'^fc'^] 

The Christoffel symbols needed to obtain Da are: 


numerical computations. The null-condition is: 

fc“fca = k°'k^ga/3 = —C^(l + 2?/')(fc*)^ + 

(1 - 2'ilj)a‘^[{k^f + r'^ik^f + sin2(6»)(fc‘^)^] = 0 

(A23) 

The four partial derivatives are: 

—2(^il)^t{k^Y ~ 2c^(l -I- 2ip)k*k*^ + 2[(1 — 2tp)aa^t— 

V',ta2][(fc’')2 ^ ^ ^2 sin2(0)(fc‘^)2] + 

2(1 — 2ip)a‘^[k^k^^ + r'^k^k^^ + sin^(0)fc‘^fc^] = 0 

(A24) 

-2cV.r(fc‘)^ - 20^(1 -h 2V')fc*fc‘^ - 2il}^ra\¥ f+ 
r'^{k^f + sin2(6»)(A:'^)2] -f 2(1 - 2'il;)a^[k^k''.^ + r{k^f+ 

r^k^k^j, + r sit?( 9){k‘^Y + sin^(0)fc'^A:'^] = 0 

’(A25) 


r‘ = 

“ 1 + 2'? 
pt _ 

1 + 2 ? 

pt _ ?,B 

1 + 2 ? 

pt ^ 

1 + 2 ? 

rr -rS _ r<^ - Q^.t + Q,t(2V’- 1) 
a(2V>-l) 

pr ^ 

2V^ - 1 

pS ^ p0 ^ ?,rr + 2?-l 
r{2? - 1 ) 

-pr p0 ^4^,B 

^ ^ ?,e sin(6l) -F {2? - 1) cos(6l) 

~ {2? - 1) sin(0) 

pr p0 p<?:> 4^+ 

^ <pr — ^ <pe — ^ — 2? — 1 


(A21) 


Using the Christoffel symbols, the following differential 
equation for the angular diameter distance is obtained: 


4V' + 1 , o“.‘7t 




u4> 
'^4 

cos( 6 ») g 


(A22) 


2? ' ; + 3^fc* + -fc" + 

4:?^ — 1 a T sin( 0 ) 


When setting the initial conditions, the null-condition 
and its partial derivatives are used to ensure that 
the geodesic will be null, the null-condition and its 
derivatives are also used for checking the accuracy of the 


—2<??fi{k*Y — 2<?{1 + 2?)k^k*0 — 2?^ea^[{k'^)'^+ 
r‘^{k4'f + slT?(9){k'^f] + 2(1 - 2?)a^{k^k'^g + r^k^k^g+ 

sin^(0)fc‘^fc g -I- cos(9) sm(9)(k'^)^) = 0 

(A26) 

-2c^?^^{k*Y - 2c^(l -I- 2?)k*k^^^ - 2?^^a?[{k'^Y + 
r‘^{k^Y + sin^(0)(fc'^)^] -|- 2(1 — 2?)a?[k''k'^,f,+ 

r'^k^k^^ + sin^(6*)fc‘^A;‘^] = 0 

(A27) 


For an observer placed at the origin, the initial 
conditions must be in accordance with a radial null- 
geodesic. Aside from the trivial constraints this sets 
on the initial conditions, two constraints are worth 
mentioning. First, the initial condition of k\ should be 
determined from the following expression: 


k 


t 



1 

\T, R, 


?,t a^t 

2 

T + A J 


1 + 2? 1 — 2? a 


(A28) 


The initial condition of k? is determined from 
the partial r-derivative of the null-condition 

with k^?o = -5 -I-^ |o = 0 , i.e. 

Kr\o = (i-tp)a'^kr + a^{k^f). 

See e.g. [25] for details on how these initial condi¬ 
tions can be derived. 
























Appendix B: ODEs for the quasi-spherical Szekeres 
model in spherical coordinates 
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In this appendix, the ODEs used to obtain Da{z) in 
quasi-spherical Szekeres models are given in expanded 
form (in spherical coordinates). 


The ODEs are solved by using a gsl ODE solver. 
In order to do this, it is necessary to eliminate and 
k‘^ from equation (24) and -^{k^a) and ^(fcJ^j) from 
equation (28). This yields two new equations which are 
given below. 

The first equation is obtained by using equations (25) 
and (26) to eliminate fc® and k'^ in equation (24). The 
resulting equation is: 


fc’’ + ^ ^ [Rk^ + 

^ P F 

-Pk^ --^k^ +—R4kn^+ 
- %Fk^- 

r r r 

r r 

^{RAkl^ + FAk^f+ 

PAFf + -k 2$,,.fc'^fc’')] = 0 


(BI) 


|0) + F(i?,„ - - |0,„) + ¥{R,^pP+ 

RA%) + 

+ e,Ai){k^ - 

’ r 

^PA + k%A - ^P) + - |f„)+ 

a(0 - 1^) - dP,rAkl^ + Rr¥k''^ + ^F™(fc^)2 + 
F,rk^k% + + P.rPPc. + Q,rA^¥ + 

eiAk%k^ + k^kA + + <^Ak%k''+ 

k^KJ) - ^kAp,cA^ + PAi) + ^i\R,Akl^+ 
RA'Ko. + + QAk%k^ + k^kA + 

^Akik^ + PkA)) - ^k\F^^pP+ 

PAi) + ti^RMAf + RA^k^a + lPMkA+ 

' r Z ’2 

PfiPk% + + e,eik%¥ + AkA+ 

+ <^>Akt¥ + k^A)) 

(B2) 


Equivalently, equations (29) and (30) are used to elimi- Inserting a = t,rAA into equations (27), (29), (30) and 
nate ^{kA and -^{kA from equation (28) to obtain: the equation above, gives 16 ODEs for fc(^. The four 
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equations for above are (expanding the /3-sums): 


+ k^,e,+r.ko^, + r^k^,)- 


dX 


1 . $ . 0 . 


D ^ ^ 

^ P F 


F 


k^iRt - - %e,t) + k^{R,ttk* + R,trk^ + R,t9k\ 

r t 

i^Mk^ + + ^,tek^ + ^Mk^ + + 

+ <i>,sk^t + ^,4,k*){k^ - + (0,„fc‘+ 

+ Q^tek^ + 0,t^fc^ + 0,iA:‘,+ 

0,,fc’, + 0,eA:^, + Q,^k\){k^ - %¥) + 

|pt) + fc;^,(<i> - |p) + fc^(0.t - |pt)+ 
(0 - 1^) - (^^,rt(fc’')" + + ^Prt(fc")" + 

+ ]^P^rt{k*f + Prfc'^A:;^ + Q,rtk^¥ + 
QAk^,tk" + + ^.rtfc^fc" + $.r(fctfc’' + fc-^feyt))- 


^kAP^ttk* + P,trk^ + P,tek^ + P,tk't + P,rAt + P,ek'^t)+ 


$ .1 


+ R,A"Kt + fc" + ^AKtk''+ 
A At) + + ^Ak^k^P 

k^kA)-^kAFttk* + Ftrk^+ 
’ r 

FA*t + F,rk:t) + ^ilRMkA+ 

R^k^At + ^PM{k^)^ + Pfik^k^t + Q^teAk-^P 

e.gikA"' P Ak\) + ^,0tAk'^ + ^AkA"' P k'*’k\))) 

(B3) 


A^At) = AAtk'r P ArKr + Aek^,r + AA^)- 


1 . $ . 0 . 

—^(F,(p--$--0)+ 


7? _ £3 _ ^ 
n p F 


p 


k''{R,r - - 50 ,r) + kAR,rtk^ + R,rrk^P 

' P ^ F ' 


R rok^ + R pfjyk^ “1“ R fk^,^ + R pk^^ H- R ^k^^F 
^,4>^tr) {^,rtk^ F ^^rpk^ + ^ ^gk^ + ^ p^pk^F 

+ <S>,rA^ + + <s>,At-)F - + 


+ Q^rpk^ F ^,r(pk^ F + 


0,,fc; + 0.eA:t + Q,At)F - 1^") + 

|P,) + A:^,(<j> - |P) + fc^(0,, - |f,) + 

fc®,(0 - |f) - {\R,rr{kA'^ P RA"Kr + + 

F,rAAr P \P,rr{kA‘^ + P,tAAt + P>,rrAk^ + 0,^(A:^^F + 

+ ^FrAk^ + + A At)) - ^ (T^.rt fc* + 

P^fc" + P.^fc® + PtA:‘, + P,fc(, + PfiAr)P 
^{\R,AkA P R,A"Kr P P>,rA'k'^ P &AArk'' + 
AkA p ^,rA^k'' + ^AArk'^ p AkA) - %kAFFtkA 
P^fc" + FA\ P F,rkA p ^ilRMkA P R,ek"ArP 

' ' r Z ' 

\pMkF P PfiAAr P Q,reAk'^ + 

0,e(fctF + AkA p ^firAk'^ p ^AArk'" + k^kA)) 

(B4) 
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jL 

dX 




-(KtK 


+ KrK 


+ KeK 


+ k: 



n p F 

k’'iR^ - - fe,^) + k’'iR^tk* + R,^rk^+ 

r r 

R,(l>ek^ + R,(j>(j>k'^ + R^tk^cj, + R^k^tp + R,dk^pd- 
R,<f>k^p) + ,(ptk^ + ^ ,<t>rkX + + 

(fc"^ — pfc*^) + {Q,<j)tk* + ^,4>rkX + 0_0efc^+ 




(fc:.) 


n _ ^ ~ ]5^~ 

R p F 

0 • • $ 0 
_0) + fc-(i^,__$,__0,)+ 


k^ (^Rgtk^ + Rfirk'^ + R.offk^ + Rfipk^ + ki^tk^fi~k 


R^r^g + 

+ ^,rKe+ 
<b,ek% + <^,pk’^g){k^ -^kn+ 


&,p>p>k* + &,tk*p + 0,.r^ + 0,efc'^ + 0,0fc'^^) 
(fc«-|F) + fc^<i>.^ + fc^^($-|p)+ 

fc^0.^ + A:®^(0 - |/^) - + 

+ Frk^k% + P,rk*k'^p + 0,,0/c®fc"+ 

0.,(fc«^r + + fc-^fc;-^))- 

+ 0,^^/c®F + e,4k%k^ + k^Kp)+ 
+ ^Ak^k^ + k^A) - jk\FA^ + F^A 

+ ^{\RM{kl^+Rfik^Kp+ 

Pfik'^k'^p, + e.^efc^fc" + 0,e(fc®0fc"+ 

fc«F^) + $.e0fc^F + + fc^fc"^))) 

(B6) 


(0,etfc* + 0,e.F + 0,eefc® + 0,90^'^+ 

0,tfc*e + QA"e + + 0.04)(fc® - %k'^)+ 

- ^PA + + fc'0,e+ 

fcV© - 1^) - {\R,re{k"f + + FFAk%+ 

+ PA'^k* + 0.^efc^F + Opik^gk^ + AkA+ 
+ k’^^s)) - ‘^k^iP,etA + PfirkA 


The four equations corresponding to equation (27) are: 

^{kA = -Atk], + kA:t + k^gAp + kA^A 
^[\RAkl‘^ + R.tk^ At + \F,u{k%‘' + F,tk'k^,t+ 
\pAkA^ + P^tAk’^t + ^^k'k^ + $,t(fc;^fc"+ 

+ QptAk'^ + 0,i(fc® F + fc^/c’))] 

’ (B7) 


Pfigk^ + P,tAe F P,rAe F PfiAo) F ~p{';^R,(t>0{kX)^F 
RA'kfi F F QAkfik'' + AkA F 

<^Ak‘!’ek^ F k^A) - ^k^FpAg + F.F,) + 

+RA"Ke F \pMkA^ F PA^Ae + ^MkA'^+ 

OAk^k'- + AkA F <^>,09^^ + $,0(4^" + AkA)) 

(B5) 


— AA = -lAtk'r F ArKr + k^k^r + k^Ar]- 

k[\R,tr{k''? F R,tk"k''r F l,FAky F F^^k^ 
c Z ' Z 

\p,tr{k*? F PpkA^^ + + $,t(fc;t^'-+ 

k^kA F e,trk'^k^ + 0,t(fcU'' + kAA] 

’ (B8) 
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^ilRMkl^ + Rtk^Ko + Ftk^k>^e+ (gg^ 
IPMk^f + P,tk^k*, + + ‘^,t{k^ek"+ 

k^¥s) + Q^tek^k^ + 0.t{fc® F + k^k^e)] 

^ik^^) = -[k*,tk*4> + KrK<P + ^'ek% + k]^k^^]- 

+ R,tk^K^ + F,tk^k<^A+ (BIO) 
P,tk*k^^ + + ^,t{k^^k'^+ 

k^¥^) + + 0.t(fc«^F + k^¥^)] 


^(fc^6») — + k^rk^^e + k^,ek%+ 

f^'A) - ^[k''iF.tk^,e + F,rK,) + Ffc^,+ 
F(0,tefc‘ + 0,.efc’' + 0,eefc® + 0.0efc‘^ + 0,t/c‘e+ 
0,,^^ + 0,0fc®g + 0,0fc^g) + ©Fg + 0,eF+ 

0( + Ktk^^e + + F.fc®, + 

[^i?,e,(F)2 + Rek^k^g + ^Pee(fc^)2 + P,ek‘^k'^g+ 
Qfiek^k^ + 0,e(fc®eF + + $,eefc‘^fc’'+ 

(B13) 


Equation (29) has been used to obtain -^(k^a)- Inserting 
a = t,r,9,(j) the following four equations for ^(k^a) 
obtained: 

^(^ft) = ~ik^tk*t + k^rKt + k^ek^t + k^,(j>k’^t)~ 

4[Ftfc® + fc®(E„fc* + F^trk^ + Ftk\ + F^rk\)+ 

Ffc® + F(0,ttfc‘ + 0,rtr + 0,,ifc« + e,^tk^+ 

0 t/c‘( + 0 rk^t + 0 efc® + 0 ^ki) + &k^t+ 

• (Bll) 

Q,tk^ + ^i^iKt) + Ktklt + KrKt + fc:,fc« + 

K^k^t) - ilPMkl^ + Rek'-Vt + \pMk^?+ 
P,ek'^k’l‘t + Qfitk^k'^ + 0,e(fc® fc'’ + k^Kt)+ 
^Mk'^k'^ + + fc-^fc;;)]] 


~^^k%) — —{k^tk*^ + + k^,ek%+ 

k%k'^4,) - lik^F^tk*^ + E.fcC^) + Fk%+ 
r (0,t^A:‘ + e,r^k^ + e^g^k^ + 0.^^/c^ + 0,tfc;^+ 
0,,F^ + Q,ek% + + 0fcC^ + 0,^^"+ 

©(^(fc:^)+ + fcc.fcc^++F^fc^^)- 

[^^^0(^)2 + Rgk^k'^^ + p,ek*k^^ + 0,e0fc^F+ 

0,e(fc^^F + fc^fc)^) + + fc^F^)]] 

(B14) 

Isolating ^(fc;^„) in equation (30) and inserting a = 
t, r, 4>, 9 leads to the following four equations: 


A(fcej = -(fc^^fc*, + k^rKr + k^ek^r + k%k^,)- 

4[F,fc« + A:®(Ertfc‘ + + FtfcV+ 

+ Ffc® + k^ie,trk* + e,rrk^ + 0,e^fcV 
0.0rfc'^ + 0.tfc*. + 0,rfc(; + &,0k^r + 0.0^:t) + 
0fc; + 0,,fc" + 0(:^(/c;) + A:'*fcV+ 

k^,k^^ + k^gk'^, + fc;^fc;t) - [iRMk'f + Rek^k^r+ 
^PMk'^f + Pfik^k’^^ + Q^srk^k'' + e,g{k^,k^ + fc®A:(;)+ 

<^.grk^k^ + $.e(fctfc’’ + /d^/c;)]] 

(B12) 


^(fc^) = -{k% + + fc;^^A:^)- 

+ fc^(P„/c‘ + PrtA:" + P,etk^ + P^^tk^ + Pt/c‘,+ 
P,rk^t + P^gk^t + P^^kF) + Pk'^t + + 

^ ,trk^ + ^^tek^ + ^ ,t4>k^ + ^,tk*^t + 4*,r-fc)t + 

+ $(-^(F,) + 

’ ’ ’ dX ' 

k^.e, + + fcC.fc® + fc;'^A:^)- 

^t(fc")2 + R^^k-^r^ + 0,0tfc®r + Q^k^k-^P 

k<^k\) + 
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+ kp^,)- 

+ fc'^(Pt,fc* + P,,F + Perfc® + P^rk'^+ 

Pp^ + P,rP + P,Bp + P^p) + Pkp 

k^i^Atk* + ^.r-rfe" + ^,rsk^ + ^ A^k'^ + ^ ,tkp 

+ $.,F+ 

) + k\k*^ + fc;fc; + k^gkp 
dX ’ ’ ’ ’ ’ ’ ’ 

KPt-) - llPMkp + R^k^Kr + + 0,0(fc® F+ 

fc^fc;) + + $.0(fctfc’' + fc’^fc;)]] 

(B16) 


kpfi) - ],[Pfik^ + + P.rsk'^ + ^.eefc® + 

P,<pek'^ + P,tk^,s P P,rP\6 P P,6k^0 P P,4'k‘^g)P 

Pk*, + fc’’($,etfc* + + $,e0fc^+ 

+ ^,ekO, + + 

<^.9k'' + + k\k\ + 

aA ’ ’ ’ ’ ’ ’ ’ 

Kp‘^g) - [\RMkP + R,<^k^K9 + Q,^ek^k^+ 

QPkp^ + fc«F,) + <^.g^k^k^+ 

<^Akp'^ + k^K9)]] 

(B17) 


— ~ikp^,,ji + prKfj, + kp^^ + kpp- 
^[PP'^ + k^{P,tp + P^r^k'' + Pfip + P^^k'^+ 
PP^ + P,rk:^ + P^gk'^^ + P,^kp + Pkp 

+ <^,tkp 

+ <i>,gk% + <^^^kp + <i>F^ + $_^F + 
$(^(fc:^) + Ktp + KrP + Ksp + Kpp- 
[\RMkP + RA"P + 

k^P) + + k^Kp] 

(B18) 
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$ 0 
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p .F 
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r® = — 

re 2P 
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$ $ 0 
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(-^.F — — 55 ^{R,(p~ 

n p p 

0 0 $ 

p0.0+p$.e-pPr)) 


The initial conditions are set using considerations 
analogous to those made in section V B of [25]. The 
result is that k*^ = —0.5(lninitially, while the 
remaining fc'^ are zero or determined from the partial 
derivatives of the null-condition. 


The Christoffel symbols needed in the ODE for 
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